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TRANSONIC AIRFOIL FLOWFIELD ANALYSIS 
USING CARTESIAN COORDINATES 

By Leland A. Carlson 
Texas A&M University 

SUMMARY 

A numerical technique for analyzing transonic airfoils is presented. The 
method employs the basic features of Jameson's iterative solution for the full 
potential equation, except that Cartesian coordinates are used rather than a 
grid which fits the airfoil, such as the conformal circle-plane or "sheared 
parabolic" coordinates that have been used previously. Comparison with pre- 
vious results shows that it is not necessary to match the computational grid to 
the airfoil surface and that accurate results can be obtained with a Cartesian 
grid for lifting supercritical airfoils. 

INTRODUCTION 

The first practical numerical technique for solving two-dimensional steady 
transonic flows was developed by Murman and Cole ^ for the small perturbation 
equations in Cartesian coordinates. The essential features of this technique 
were the utilization of retarded (upwind) differences in regions of supersonic 
flow, central differences in subsonic zones, and the application of a far field 
boundary condition at a finite distance from the airfoil. The surface flow- 
tangency boundary condition was applied on the horizontal axis, consistent with 
the small disturbance approximation; and the resultant finite difference equa- 
tions were solved iteratively by numerical relaxation. Shock waves occurred 
naturally in the course of solution as steep compressions smeared over several 
mesh points, and accurate results were obtained for slender airfoils. 

This technique was subsequently extended to the complete potential flow 

(?) 

equations and the exact boundary conditions by Garabedian and Korn v , Steger 
and Lomax ^) , and Jameson ^ and was successfully applied to thick, blunt- 
nosed, and aft-cambered airfoils. In these works the edge of the computational 
coordinate system was closely matched and aligned with the airfoil surface in 
order to accurately represent the surface flow-tangency boundary condition. 

For example, Steger and Lomax used a curvilinear wraparound system, while 



Garabedian and Korn and also Jameson conformally mapped the exterior of the 

airfoil onto the interior of a unit circle with infinity corresponding to the 

(41 

circle origin. Jameson also used a parabolic type system ' ' in another 
version designed to handle supersonic freestreams. These transformations had 
the advantage that near the body they aligned the coordinates with the flow, 
which is particularly advantageous in the leading edge region where body slopes 
and vertical velocities can be large. Unfortunately, such transformations 
introduce into the differential equation many, sometimes complicated, trans- 
formation derivatives. 

Now in order for any backward difference scheme at supersonic points to 
work properly, the coordinate along which the backward differences are taken 
must be closely aligned with the flow direction. This alignment is very 
difficult to achieve with mapped coordinates when the freestream Mach number 
is high or the supersonic zone is large. On the other hand, improper alignment 
at supersonic points can lead to situations where the velocity component in 

each coordinate direction is subsonic, that is 

22222 

U,V <a <U +V (1) 

This situation not only can lead to numerical instability but also can impart 
to the finite difference equations an incorrect zone of dependence. To remedy 
this problem South and Jameson have introduced a rotated finite differ- 

ence scheme which simulates a local rotation to coordinates along and normal 
to the velocity vector; and, in so doing, they have created a scheme that not 
only has the correct zone of dependence but also does not require coordinate 
alignment with the flow. 

As a consequence of this development, the questions arose— Can the full 
potential flow equations be solved for a transonic case in a simple Cartesian 
or stretched Cartesian system by using the concept of rotated differences? If 
so, what is the accuracy and what are the properties of such an approach? In 
order to answer these questions, the present study was initiated. Its primary 
objectives have been 

(a) To develop an analysis program using Cartesian coordinates 
(airfoil given, Cp and flowfield unknown) 
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(b) To develop a design program in Cartesian coordinates 
(Cp given, airfoil and flowfield unknown) 

(c) To study the results obtained by such an approach and 
to determine the accuracy and idiosyncrasies associated 
with a Cartesian system. 

This report describes the development of, and results obtained with, the 
analysis program. Reference 6 discusses the design program. 
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SYMBOLS 

Coordinate stretching constants 
isentropic speed of sound 
coordinate stretching constants 
pressure coefficient, (p-p m )/(h pJJ^ 2 ) 
lift coefficient 

coefficient of moment about the leading edge 
coordinate stretching functions 
Mach number 

coordinate normal to the streamwise direction 

pressure 

vel oci ty 

streamwise coordinate 
time 

velocity component in the x-,y- direction respectively 
Cartesian coordinates 

angle of attack or coefficient in Eq. (34) 

1-M 2 or coefficient in Eq. (34) 

ratio of specific heats or coefficient in Eq. (34) 

circulation 

damping coefficient 

polar coordinate 

computational coordinates 

densi ty 

potential function, Eq. (3) 
perturbation potential, Eq. (4) 
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Subscripts : 


b,body 
TE 
i J 

€,n,x,y 
Superscripts : 


freestream condition 
body 

trailing edge 
grid location 
differentiation, i.e., f 


af 

ax 


+ Value computed by current relaxation sweep 

PROBLEM FORMULATION 

The exact equation for the potential function for two dimensional 
pressible flow can be written in Cartesian coordinates as 


.2 2 2 2 

(a - $ )$ - 2$ $ $ + (a - $ )$ =0 

' x 1 xx x y xy v y ' yy 


where 


_ a$ _ d*± . dh 

x ax, xx gx 2 , $ xy = axay 


etc. 


If a perturbation potential, <p, is introduced such that 

o = xq m cosa + yq^sincx + qj> 
where the velocity components are given by 

U = $ x = qjcosa + <i> x ) 

V = $ y = qJ sina + ty) 

then the governing equation for the perturbation potential becomes 


wi th 


(a 2 - u\ x - 2UV* xy + (a 2 - V 2 )* yy = 0 


2 2 , v _, 2 2 2 - 
a = a - ( V- > ® * v - 3 


The pressure coefficient at any point is given by 


r = P ~ P-” = —i—o f n + X~J- m ( 1 
u p u 2 ^r 2 1 L1 2 » u 

' ^ CO CD 1 CO 


„ 2 2 

2 - u +v 


Y-l 

2 -)] - 1 } 


com- 

(?) 

(3) 

(4a) 

(4b) 

(5) 

( 6 ) 

(7) 
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The appropriate boundary conditions are, at the airfoil surface 


(^■) = ( tt ) 

ax body u body (8) 

and at infinity (2»7,8) 

<P = Tan -1 (3 Tan(e - a)) (9) 

where 8 is the polar angle and r is the circulation determined by the change in 
potential across the Kutta-Joukowski cut at the trailing edge of the airfoil, 
i.e. 


r ^y=0 + " ^y=0“^Trailing Edge 


( 10 ) 


NUMERICAL ANALYSIS 


Coordinate Stretching 

In the present problem an infinity boundary condition must be applied at 
the edge of the computational grid. It is convenient to transform or stretch 
the original x,y plane to some 5,n plane so that the edges of the ?,n grid 
correspond to infinity. In this way, Equation (9) can be applied directly. If 
this is not done, a far-field condition can be used at the grid edges 
instead. 

After investigating several possibilities, the coordinate stretching 
represented on Figure 1 was selected. Here x,y is the physical plane and £,n 
represents the computational plane, and each is subdivided into three regions. 
The stretching is symmetrical about the origin and it is given by 

x = x 4 + A 2 Tan[ \ ($ - S 4 )] + A 3 Tan[ yU - 5 4 ) 3 ] ( U ) 

in region III and by 

2 

x = §(a + b 5 ) ( 12 ) 

in region II. The constants a and b are determined by the requirements 


and 


x = x 4 


at C = £ 

4 


(13a) 


dx _ ttA 

d£ ~ 2 


at 5=S 4 


(13b) 
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The constant controls the grid spacing in the vicinity of x^, usually near 
the leading and trailing edges; while determines the physical location of 
the grid line adjacent to the grid edge. 

In the y-direction the stretching relationship is given by 

y = A x Tan( \ n) (14) 

where Aj controls the grid size near the airfoil. 

Notice that the stretchings map the infinite physical x,y plane 

-00 ^ X - 00 

—oo ^ y 5 co 

into the finite computational plane 

-(1 + 5 4 ) S 5 < 1 + s 4 . 

-1 < n - 1 (15) 

where determines the amount of the computational plane confined to the 
vicinity of the airfoil. A typical grid system is shown on Figure 2. 

When selecting a stretching care must be taken to insure that the 
stretching does not force a physically unrealistic or abnormal behavior on the 
solution. For example, analysis of a compressible doublet indicates that <j> y 
and 4> v should decay as 

A 


-y' 3 

as 

y -»■ 00 


- X - 3 

as 

X -> 00 


- X - 2 

as 

x » 


_ -2 
y 

as 

y -* 00 

(16) 


In the present £,n system, this type of decay would require that 

<i> ~ y" 1 as y ■+ » 

~ constant as x -*■ °° (17) 

and this type of behavior is quite permissible. However, with some other 
stretchings such as n=tanhy, the decay indicated by Equations (16) would re- 
quire that 

<j> ~ « as y ->■■«>, n 1 (18) 
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In practice, it usually happens that the derivatives in the computational 
plane remain finite, even with the exponential stretchings; but this precludes 
the correct asympototic decay, and the effect may produce stable, convergent, 
wrong solutions. 

Now as a result of the introduction of coordinate stretching, the 
governing equations should be written in terms of the independent computational 
variables £ and n. By defining 


f = 

Equations (4) and (5) become 


d£ 

dx. 



(a 2 - U 2 )f(f<j> c ) 5 - 2UVfg^ n + ( a 2 - V 2 )g(g* n ) n = 0 


(19) 

(20) 


U = qjcosct + f* 5 ) 

V = qjsina + g<f> n ) (2 i) 

Finite Difference Scheme 

As indicated in the Introduction, in order to avoid at supersonic points 

difficulties associated with nonalignment of the coordinates and the flowfield 

a rotated finite difference scheme is used in the present problem. In this 
(4 5l 

approach v ’ ' the principal part of the governing differential equation can 
be written in coordinates parallel and perpendicular to the local velocity 
vector, S and N respectively, as 

(1 ' h ) Hs + % = 0 (22) 

d 

where 

4>ss = U 2 f(f^) ? + 2UVfg* ?n + V^g*^] 

4> NN = V 2 f(f^) 5 - 2UVfg^ ?n + U 2 g(g<}» n ) n ] .(23) 

Notice that substitution of Equations (23) into (22) yields the governing 
equation, Eq. (20); and thus, as pointed out by South, ^ Eq. (22) is simply 
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a rearrangement of terms which exhibits the basic features of a local rotation 

to the streamline direction. 

(4 5) 

The basic concept ' * ' of the rotated scheme is to use, at supersonic 
points, first-order, upwind differencing in both the ? and n directions for all 
contributions to ^ and central differencing for all contributions to <j>^. In 
this manner the correct zone of dependence at supersonic points is built into 
the finite difference scheme. At subsonic points, the normal procedure of using 
central differences for all derivatives directly in Eq. (20) is used. Conse- 
quently, the scheme is second order accurate at subsonic points and first order 
accurate at supersonic points. The resulting finite difference equations are 
solved iteratively by using column relaxation sweeping from upstream to down- 
stream. In such a procedure the values obtained during an iterative sweep can 
be thought of as new values, while those obtained during the previous sweep 

• J 

can be considered old values, $.... In this manner, the change from one itera- 
tion to the next can be viewed as that occurring during some artificial time 
step. At; and artificial time derivatives such as <f> t , ^ and ^ can be consid- 
ered. Of course, as the relaxation process converges, these terms become 
negligible. 

Various investigators have used different finite difference formulas to re- 
present the derivatives in Eq. (23). For example, Jameson ^ and South and 
(51 

Jameson ' ' used in their finite difference formulation a mixture of old and new 

values. In this way, terms such as and <j>^ were implicitly added to the 

problem and used to control numerical stability. Since it was determined by 
f4l 

Jameson ' ' that additional <j>^ t needed to be explicitly added to the problem in 
order to insure stability, it was decided in the present case to use all old 
values in the difference expressions for ^ and to add explicitly. For 
appropriate combinations of old and new values were used, causing <j>^ t to 
occur implicitly. It should be emphasized again. that <t> Nt and <t> st go to zero 
as the solution converges, and hence do not affect the final result. 

Thus the finite difference expressions used in the present formulation are: 

2 2 

For contributions to <(> MM when q > a 


(1 V? = ^‘ { W*i+i,j • ‘NV • W4j ■ 4-i, j )} 


(24a) 


4A£Ari 


. T T 

^i-i.j " +i-ij+i “ ♦i+i.j-i < h'+i,j+i} 
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(24c) 


(g<fc ) = {g ... (<j>t ... - - g . . (<j)t . - *1" . , )} 

Vn .2 ■ m j+V v i,j+i U s J-V v iJ 


An 

2 2 

For contributions to <p<.$ when q > a and V > 0 


{f Vc = a ^2 {f i-% (<f, ij ' +i-i,j ) " f i- 3 / 2 ( VlJ “ *i- 2 ,j )} 


1 


Co A?An 


(A. . - d> . . — d> . . +<b. . } 

1J 1-1, J 


^n^n Ati 2 " g j-3/2^i,j-i " 

where to create stability there is added to the basic equation 


-eAt f _ -efAt rUf . Vg 

~ AT *€t "q +nt ] 


-ef Uf i -1 
AC ^ A?q' ' Y i j T i j 


((f).. “(f).. — (j). . ■f’d). .) 


+ y. ... izk . (*t. - <t. . - *1 . + <$>. . )] 


For contributions to 4 > ss when q > a and V < 0 
(1 V§ = ‘ ♦i-U* " f i-3/2^i-i ,j ‘ *i-2,j )} 

<t> = — { -d> . . + <fc . . + 6 . - 6 . } 

v Cn ACAn 9 iJ m-i ,J v i,J+i y i-i,j+r 


^n^n A 2 ^j+J^ij " ^ij+i^ " g j+3/2^i J+i “ <f> i,j+2^ } 
An 

and to create stability there is added 


-efAt 


= -if [U _Aih - * + . + .) 

L q A? VH ij MJ M-i,J M-i.J 


AC y St Ac L q AC 

+ i h+h . + $ + - 4 ,. )] 

q An v 9 lj MJ M,J+i 


(25a) 

(25b) 

(25c) 


(25d) 

(26a) 

(25b) 

(26c) 


(26d) 
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2 2 

At points where q < a the expressions are 

+ 

(f Vc = (f Wi+i,j ' (f i+H + f i C “W 1 + 

+ 

_ 1 . + + 

*5n " 4A?an {<f> i-i,j-i " *i-i,j+i " *i+i ,j-i ' *i+i ,j+i} 


(27a) 

(27b) 


(gd> ) = — {9 -j_i (<f>* -$*•)- 9- > ( - - d>t . )} 

n n A 2 y J+% V M,J +1 9 lj' a J-VMJ 


(27c) 


where the relaxation factor, w, has been incorporated into the difference 
formulas. 

In all cases the velocities, U and V, are represented by central diffe- 
rences using old values. When these expressions, Eqs. (24) to (27), are sub- 
stituted appropriately into Eqs. (22) and (23) or Eq. (20), the result is a 
tridiagonal system of equations that can be solved for the values of <j>^ • on 
col umn i . 


Treatment of Boundary Conditions 

There are many ways to approximate the flow tangency condition, Eq. (8), 
at the airfoil boundary. One of the simplest, which is used here, is to 
generate dummy values of <f> at mesh points inside the boundary such that the 
usual difference equations can be solved at points just outside the boundary. 
The problem is to develop a scheme for providing and updating these dummy 
values by using the surface flow tangency condition and neighboring values of 
<j> in the mesh, with adequate accuracy and without creating instability. 

To accomplish this, we first note that, in the computational coordinates, 
Eq. (8) becomes: 


/dyx _ ^b _ sing + 9 b^nb 

'dx U. , f. <j>_. 

b b cosa + b £b 


(28) 


Next a Taylor series about the dummy point (i,j-l) is written (see figure 3): 


4 , = K 
C b 5 i ,j-i 


+ (n b • j., + -‘- 

(29a) 


(29b) 
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When these are written in finite difference form using second order expressions 

for 4> and <f> and at least first order ones for <j> and <j>,. they become (for 
5 n T on So J v 

the upper surface case) 



-3d> . . +4 d> . . - <f> . . 

1,J-1 y lJ »J+l 

2&0 



V)>< 


♦i-M ' 


An 


(30a) 




i+U-i ~ 

2A 5 


+( n - n )( t l ± Lj ~ iliLJil ) 

n j-l M 2 ASAo ; 


(30b) 


These expressions can then be substituted into Eq. (28) and the result solved 

for a sufficiently accurate . that is in terms of the neighboring po- 

1 » J " * 

tentials, body slope, and body position. 

For those situations where the flow at point (i,j) is supersonic and 
V.jj > 0 examination of Eqs. (25) show that the difference scheme will require 
a <b. . value as well as <j>. • Thus, in all cases a value of <j>. . is 
determined by extrapolation as 




i ,J~z 


-6. . + 2d>. • 


(31) 


The above procedure for determining the values of the dummy mesh points 
inside the boundary is performed twice for the relaxation procedure of column 
i. First, Eq. (28) is used with ()> values obtained by the previous relaxation 
sweep in order to obtain old values for the dummy points 4 • . and <j>. . . 

1 9 J~Z 

Then, after column i has been relaxed, it is used again with as many current 

values of <j> as possible to obtain new values <|>t . and <j>* . . In this 

* j J 1 2 

manner, the dummy mesh points will have both old and new values just like re- 
gular mesh points, and they can be used directly in the finite difference 
formulas without special treatment. 

A similar procedure is used to satisfy the boundary conditions on the 
lower surface of the airfoil. 

At the edge of the grid the infinity boundary condition given by Eq. (9) 
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must be used. In the stretched £,n system the polar angle will have at the 
edge of the grid the values shown on Figure 4. Thus, the boundary values will 


be 

On line AB, 0 = 0, <)> = (-Tan -1 (aTana)) (32a) 

On line BC, e = j, <j> = Tan” 1 (ecota) (32b) 

On line CD, 6 = it, <p = ^ (tt - Tan" 1 (eTana)) (32c) 

On line DE, e = J^, = §“ (* + Tan" 1 (ecota)) (32d) 

On line EF, e = 2tt, <j> = (2tt - Tan” 1 (eTana)) (32e) 


Boundary values must also be assigned to the singular corner points B, C, D, 
and E since they will be used in the finite difference formula for <f>^. The 
choice of these values is somewhat arbitrary since they depend upon how the 
corner points are approached. Fortunately, numerical studies indicate that the 
solution is insensitive to the values selected; and, thus, the following have 
been used. 


Point B, * = g Tan’ 1 ( ^ ) 

Point C, ♦ -§£<» - Tan' 1 ( ^ $) 

Point D, * = g- (t + Tan" 1 ( g 

A A 

Point E, <p = ^ (2tt - Tan" 1 (^ ^-)) 

n 2. n 


(33a) 

(33b) 

(33c) 

(33d) 


These assume that the corner points are approached along a diagonal (in the ?-n 
plane) connecting the interior point just inside the corner and the corner 
points . 

Using the above finite difference formulas and boundary conditions, the 
governing differential equation can be solved iteratively by column relaxation. 
For each iteration, the sweep pattern is from upstream to downstream in the 
region in front of the airfoil, then in the region above the airfoil, then 
below the airfoil, and finally in the portion downstream of the trailing edge. 
The latter region contains the Kutta-Joukowski cut, and the equations are 
appropriately differenced to give the proper jump in the potential. 
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NUMERICAL STABILITY 


Supersonic Points 

The incorporation into the finite difference formulas of both old, <p, 
and new, <p , values and the explicit inclusion of ^ at supersonic points 
introduces time-like derivatives into the equation being solved. Thus, during 
the iterative process the actual equation at supersonic points is of the form 


n t + ( !)♦ SS + 2c “»St ‘ *NN + 2e *Nt = 0 (34) 

d 

Jameson ^ has shown that for numerical stability the coefficients must be 
such that 

Y = 0 (35a) 

a 2 > s 2 - 1) (35b) 

a 2 


In the present problem, a proper value for e can be determined from Eq. (35b) 
once a and 6 are determined. Since the 4> St is included explicitly (see Eqs. 
(25d) and (26d)), the a is known, and the B can be determined from the term. 

By adding and subtracting old values of <f> to get <f> NN in terms of old 
values, the time-like derivatives can be isolated. Thus 


^NN " 2 ^ m *- + 


Nt V NN 


V 2 f r ~ f i-% ' <|, ij ) +f (<f> i - 1 . j " 

q 2 A? 2 1-Js A£ 2 


- 2UVfg r 'N-iJ-i " ‘N’-i.j-i ~ ‘N-i.j+i + ^i- 


[ 


4ACAn 


ill »J +1 J 


+ U 2 q r a ^i J+i ~ " ^i j + *ij^ 

2 9 i+Js . 2 

q L An 


($ + . 


■H 


JA 


*i.j-i + 
An 2 


.J 


J 


J 


(36) 
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( 36 ) 


= <f>, 


V f At 


NN 


At ft. v , UVfg At , x , 

A? f< ^5t i-^,j 2 AC ,j 


* U f &K" Vi j* -<«♦„* W 

which in the limit of vanishing step size becomes 


*NN - 2B *Nt * *NN + T « [ ? f *Ct + f (37) 

where the term in brackets is Hence., from a time-like viewpoint the actual 

differential equation being solved is 


(1 


4 * 


eAt 

SS ' AC 


f »st + 


+ l k 4t *Nt = 0 


and so 


Y = 0 


(38) 


- eAtf . _ -Vf At 

a 2AC S qAC (39) 

Thus, from Eq. (35), a necessary condition for stability is 



Notice that this requirement means that, for a fixed value for e, 
numerical instability is most likely to occur where the local Mach number is 
large, say immediately upstream of a shock wave. Numerical experiments have 
been conducted to verify Eq. (40), and in most cases instability, if present, 
does originate from the high Mach number region in front of a shock wave. 
However, the minimum value of e predicted by Eq. (40) is usually much smaller 
than the value actually required in practice. This latter phenomena has also 
been observed by other investigators using time-like damping, and usually the 
actual value is only slightly less than e v where 

iTlaX 


e max ^ ^ax 


(41) 
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Subsonic Damping 

Since at subsonic points Eq. (20) is used while at supersonic points Eq. 
(22) is solved, there is the possibility that during the initial iterative 
stages that the time dependent terms of the two finite difference equation 
forms do not match near the sonic line. This mismatch possibility can be 
prevented by adding to the subsonic finite difference equation a term 


- eAtf _ - eAtf r Uf . V _ , -i 
A? ^St AC *- q ^£t q 9 ^nt^ 


(42) 


which when q equals a causes the supersonic and subsonic difference equations 
to match exactly. 

The question naturally arises, however, of what effect does Eq. (42) have 
upon convergence and stability away from sonic points. By substituting Eqs. 
(27) into Eq. (20) and isolating the time-like terms, it can be determined 
that the actual equation being solved at subsonic points is 

2 


(1 - f < f *c>{ - ™ f 9+cn + (1 


?■> 9 (9 Vn 


. ,, U 2 . f at ,w-2, , ,, u\ f 2 ... 

A (_ir) *t ‘ (1 ' ST at * 


St 


+ M£ 9 *t , E atf [Uf + v 

2 AC 9 nt A? L q v ct q a %t J 


(43) 


Far away from the supersonic zone and the airfoil, V will normally be very 
small and the effective equation will be 

2 


(1 - \) f + 9 (94> n ) n + (1 - ~) f ~~T At (*=£■) 4> 


2 ' 2 W ' *t 

a AC 


_ [(i _ y_) lal + Mtf uf j ^ = o 

U a 2 J AC AC q J 9 ct 


(44) 


In order to see the effect of the term on the solution behavior of 


this equation assume 


f -ss g 
At ~ AC 
U ^ q 


(45) 
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and let 


A = 1 - 


^ constant 


(46) 


Then Eq. (44) can be modeled as 


A* + 


nn 


* if {l ir L) *t + ft *Et + E * s t 


(47) 


The convergence behavior of this equation can be studied following the 

(9) 

approach of Garabedian ' . First, transform from (?,n 5 t) to (5,n,s) by 

s » t + <^> E 

and assume solution of the form <}>= U(c,n)W(s). 

Thus, in general, the solution to Eq. (47) is 


(48) 


= U Q (?,n) +J. W m U m 


(49) 


where W m and U m are eigenfunctions and 




and 


_ -3 ± [p 


m 


1,2 

and C and are the Fourier coefficients, 
m m 


2a 


2 u 

4aKm ] 2 


(50) 


(51) 


Obviously, the rate of decay of the time terms is determined by the real 
part of Pjj, and the optimum convergence rate will occur when the relaxation 
factor, w, is such that the radical in Eq. (51) is zero. When it is zero or 
pure imaginary, the rate of convergence will be determined by 


= UL = -2A (J 


~vr 


w 


2a A? (A + e) ' 


(52) 


For a fixed relaxation factor, as the damping coefficient, e, increases, the 
radical in Eq. (51) eventually goes negative; and the error, E, due to N 
relaxation iterations will be of the order 


E = 0 ( e " pNAt ) 


(53) 
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p 

Since p decreases as (A + e) , an increase in e from zero to one can cause 
an order of magnitude increase in the number of iterations required for con- 
vergence. 

From this analysis it can be concluded that the addition of at sub- 
sonic points decreases significantly the rate of convergence and that large 
values of e should be avoided. 

In addition. South has shown by a linear stability analysis that for the 
model equation 

<j> + d> — eU<fc . - eV<t> , = 0 (54 

v xx ^yy v xt v yt 

that the inclusion of explicit time-like damping in addition to that impli- 
citly obtained by the finite difference forms of <j> and a can, in a 

xx yy 

stretched coordinate system, be destabilizing and may cause numerical insta- 
bility, particularly in regions of hard stretching where V is large. 

Based upon these two analyses, both of which have been verified numeri- 
cally, it is believed the specific addition of e<ji st at subsonic points should 
be avoided and that the mismatch in time- like terms near sonic points should 
be prevented by changing the damping coefficient at supersonic points to 
[M 2 - 1 f 2 e. 


NUMERICAL STUDIES 

In actually carrying out the numerical solution the parameters Aj, A,,, 
and A^, should be selected so that the grid spacing ax and Ay is small near 
the airfoil, and x^ should be near the trailing edge. The question is How 
small and How close? See Figure 5. 

To answer these questions a series of numerical experiments were conducted 
and it was determined that accurate results for lift and surface pressures 
could be obtained when 

2r 

(55) 

c c c 

in the vicinity of the leading edge. (In the present computer formulation the 
chord c is always unity. The airfoil leading edge is at x = -0.5 and the 
trailing edge at x = +0.5.) No significant changes in lift or surface pres- 
sures were observed for smaller step sizes as long as ax was about the same 
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magnitude as Ay near the leading edge. However, when Ay was three times ax 
near the leading edge good accuracy was not acheived until 

, Ay ^ 3AX (56) 

and considerable sensitivity to the value of Ax was observed. In addition, the 
pressure coefficient near the nose exhibited slight oscillatory behavior. 

Considerable sensitivity was also observed for the value chosen for x^, 
which essentially determines the first and last grid points on the airfoil. 

For a grid spacing given by Eq. (55), x^ values of 0.50 and 0.49 yielded 
pressure coefficient results that differed by as much as ten percent, with the 
results associated with the 0.5 value being worse. In addition, sensitivity 
to grid spacing was observed in the 0.5 case but was minimal in the 0.49 case. 
Further investigation showed that if the spacing near the leading edge were 
reduced to r/c, either by using more points or by coordinate stretching, that 
good results insensitive to the x^ choice could be obtained even for x^ values 
near 0.5. The cause of this behavior is not at present completely understood. 

Based upon these results and those previously discussed concerning 
numerical stability, it is believed that accurate solutions can be obtained 
when 

X 4 * 0.49 ^ 0.495 

Ax „ Ay. „ 2r LE 

c c c near the leading edge (57) 

In addition, considering the destabilizing effects of hard stretching, it is 
recommended that whenever possible Ax be kept relatively constant between the 
leading and trailing edges of the airfoil. 

COMPARISONS AND TYPICAL RESULTS 

Any new numerical technique can only be verified by comparing its results 
with those previously obtained by other investigators. As suggested by 
Lock, the NACA 0012 airfoil is an excellent test case because its shape 
can be exactly prescribed analytically; and the solutions of Sells, which 
are believed to have an accuracy to about \% , can be used for comparison for 
subcritical cases. 
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Figure 6 compares the results obtained by the present method with those 
of Sells for a nonlifting subcritical case. As can be seen the two pressure 
distributions agree exactly with the exception of the immediate vicinity of 
the trailing edge. The latter small disagreement may be due to the fact that 
in the Sells case the airfoil is extended 0.89% to yield zero trailing edge 
thickness while in the present case it is not. Also, in the present case the 
symmetry present in the nonlifting flow is not imposed, and the circulation is 
actually a calculated zero. 

Pressure distribution results for a lifting subcritical NACA 0012 flow 
are shown on Figure 7 and again the results are compared with those obtained 
using the Sells method. Here, the two sets of data are always within two 

percent of each other and the lift coefficients agree exactly. In particular, 
notice the excellent agreement on the magnitude and location of the upper sur- 
face pressure peak. 

The above cases are excellent for validating the present Cartesian 
coordinate solution on airfoils having large leading edge radii. However, 
many airfoils of practical interest have very small leading edge radii and 
extremly rapid upper surface expansion near the leading edge. For example, the 
NACA 63A006 airfoil has a leading edge radius only one-sixth that of the NACA 
0012. Results for this airfoil at M of 0.6 and an angle of attack of two 
degrees are shown on Figure 8. The present solution was obtained with a medium 
sized grid of 49 by 49, which yielded 66 points on the airfoil surfaces. 
Superimposed on this pressure distribution is a result obtained by Newman at 
NASA Langley using a program developed by Jameson. This Jameson program, which 
uses conformal coordinates in the circle plane, is believed to be very accurate 
at both subcritical and supercritical speeds, and the results shown on Figure 8 
were obtained with a 192 x 32 grid with 192 points on the airfoil surfaces. 
While the Cartesian grid had insufficient points to completely resolve the very 
sharp pressure peak, the Cp values at given x locations agree with those by 
Jameson at the same x locations. In addition, the concavity variations in the 
upper surface pressure coefficient distribution are accurately reproduced; and 
the computed lift and moment coefficients agree within less than two percent. 

Based upon the comparisons shown on Figures 6-8, it is believed that for 
subcritical cases the present Cartesian system method can yield accurate 
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results using medium-sized computational grids. For supercritical cases, 
comparison and verification is somewhat more difficult. However, comparison 
with results obtained by Jameson's program for an analytically defined airfoil, 
such as an NACA 0012, should suffice. 

Figure 9 shows such a comparison for an NACA 0012 at a frees tream Mach 
number of 0.75 and two degree angle of attack. Again the lift and moment 
coefficients essentially agree exactly, and the pressure coefficients and 
shock location agree quite well. The present results do show a slight 
re-expansion immediately downstream of the shock wave that becomes more 
pronounced with further reduction in grid spacing. It is believed that this 
expansion is detected because for this case the Cartesian grid is better 
aligned with the shock that the conformal grid. 

These comparisons demonstrate that accurate transonic flow solutions can 
be obtained in a Cartesian grid system. Further, they show that there is 
little loss of accuracy associated with not mapping the computational grid to 
match the airfoil surface and that excellent results can be efficiently 
obtained using only medium grids. 

In addition to the above comparison cases, the present analysis program 

has been applied to many airfoils specifically designed for transonic flight. 

For example. Figure 10 shows the pressure distribution obtained by the present 

method for an aft-cambered supercritical airfoil. This airfoil was 

designed ^ with a leading edge radius of 0.0245c and a weak shock ( upstream 

Mach number of about 1.20) near the three-quarter chord point, in accordance 

• . (121 

with the design criteria advocated by Whitcomb. ' ' The airfoil shape por- 

trayed on Figure 10 is not the actual surface but the displacement surface, 
and a preliminary analysis indicates that the boundary layer is free from 
separation effects. 

Figure 11 shows the result of a direct analysis with the present program 
of an entirely different type of supercritical airfoil. This airfoil was 
designed ^ to have the same basic lift coefficient and lower surface pressure 
distribution as that of an NACA 0012 and the moment coefficient was to be lower 
than that of an aft-cambered airfoil. As can be seen from the present direct 
analysis, these objectives were essentially achieved, although the lift coef- 
ficient is slightly higher. Again, the airfoil shape portrayed is the 
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displacement surface, and it has a (t/c) „ of 0.1265 at the 0.3753 chord point 

max 

with a leading edge radius of 0.0146. 

The results shown on Figures 6 to 11 were typically obtained using a 
49 x 49 grid with 66 points on the airfoil. In the actual computational 
process, a solution was usually first obtained on a 13 x 13 grid with a maximum 
of 800 relaxation sweeps, then on a 25 x 25 grid with a maximum of 400 cycles, 
and finally on a 49 x 49 grid with a maximum of 200 iterations. The converged 
results from the coarser grid are distributed over the next (finer) mesh by 
interpolation to start the iterations on the next mesh. This grid halving 
technique accelerates convergence and significantly reduces computational time 
as demonstrated in Reference 6. The convergence criteria was usually that the 
maximum change in perturbation potential be less than 1. x 10"^. Normally 
each problem required only about 15 minutes on the IBM 360/65 system or 5 
minutes on a CDC 6600 machine. 


CONCLUDING REMARKS 

A method suitable for analyzing the transonic flow about two-dimensional 
airfoils has been developed using a stretched Cartesian coordinate system. 

This method utilizes the full inviscid potential flow equation and exact 
boundary conditions and solves them via numerical relaxation. A rotated finite 
difference scheme is employed in order to obtain the correct domain of 
dependence in supersonic regions; and Jameson's time-like damping is included 
to ensure numerical stability. In addition, grid halving is used to achieve 
computational efficiency. 

It has been demonstrated by comparison with previous results that it is 
not necessary to match the computational grid to the airfoil surface and that 
accurate results can be obtained with the Cartesian grid. Best accuracy is 
obtained when the mesh size in the leading edge region is nearly square and 
the order of the leading-edge radius. This accuracy has been demonstrated for 
both lifting and nonlifting cases at subcritical and supercritical speeds. 

The program has also been successfully applied to aft-cambered and other air- 
foils specifically designed for transonic flight. 

Texas A&M University 
Department of Aerospace Engineering 
20 November, 1974 College Station, Texas 
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Figure 2--Typical Grid System (25x25) 
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Figure 3— Relationship Between Airfoil and Grid (Upper Surface) 
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Figure 6 — Comparison with Sells for Subcritical Nonlifting Case 
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Figure 9— Comparison with Jameson for Supercri tical Lifting Case 
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Figure ll--Pressure Distribution for a Supercritical Airfoil 
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